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I. The ACOR numerical code 

R-M. Ouazzani 1 ' 2 , M-A. Dupret 2 , D. R. Reese 2 



(N 

o 

Oh 

c/3 



(N 

6 

o3 



1 LESIA, UMR8109, Universite Pierre et Marie Curie, Universite Denis Diderot, Observatoire de Paris, 92195 Meudon 
Cedex, France 

e-mail: rhita-mar ia . ouazzani@obspm . f r 

2 Institut d'Astrophysique et de Geophysique de l'Universite de Liege, Allee du 6 Aout 17, 4000 Liege, Belgium 

Received September 15, 1996; accepted March 16, 1997 

ABSTRACT 

Context. Very high precision seismic space missions such as CoRoT and Kepler provide the means of testing the modeling 
of transport processes in stellar interiors. For some stars, such as solar-like and red giant stars, a rotational splitting 
is measured. However, in order to fully exploit these splittings and constrain the rotation profile, one needs to be able 
to calculate them accurately. For some other stars, such as 5 Scuti and Be stars, for instance, the observed pulsation 
spectra are modified by rotation to such an extent that a perturbative treatment of the effects of rotation is no longer 
valid. 

Aims. We present here a new two-dimensional non-perturbative code, called ACOR (Adiabatic Code of Oscillation 
including Rotation) which allows us to compute adiabatic non-radial pulsations of rotating stars, without making any 
assumptions on the sphericity of the star, the fluid properties (i.e. baroclinicity) or the rotation profile. 
Methods. The 2D non-perturbative calculations fully take into account the centrifugal distortion of the star and include 
the full influence of the Coriolis acceleration. The numerical method is based on a spectral approach for the angular 
part of the modes, and a fourth-order finite differences approach for the radial part. 

Results. We test and evaluate the accuracy of the calculations by comparing them with those coming from TOP ( Two- 
dimensional Oscillation Program) for the same polytropic models. We illustrate the effects of rapid rotation on stellar 
pulsations through the phenomenon of avoided crossings. 

Conclusions. As shown by the comparison with TOP for simple models, the code is stable, and gives accurate results 
up to near-critical rotation rates. 
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1. Introduction 

Rotation plays a key role in the evolution of stars across 
the Hertzsprung-Russell diagram (which shows the distri- 
bution of stars in the luminosity versus effective tempera- 
ture plane). For instance, the centrifugal distortion breaks 
the thermal equilibrium, and provokes large sc ale currents 

These cur- 
ities modify 



called meridional circulation ( Eddingtoii}|1925 
rents as well as shear and baroclinic instabi 



the angular momentum distribution ( jZahnj [1992 Mathis 

roil" 



|et al. 20041, and thereby the rotation profile inside the 
stars. Meanwhile, these processes transport chemical ele- 
ments, and change the evolution of the stars. That is the 
reason why determining rotational profiles inside stars is 
crucial for modeling stellar structure and evolution. 

Asteroseismology is currently the only tool that allows 
such a determination. But rotation also changes stellar pul- 
sations. The centrifugal force distorts the resonant cav- 
ity of the pulsations, and the Coriolis force modifies the 
modes' dynamic. Usually, for slow rotators, the effects of 
rotation are taken into account as a perturbation of pulsa- 



tions (see for instance 


Ledoux 


1951 


for first-order effects, 


|Gough & Thompson|1990|and 


Dziembowski & Goode|1992| 


for second-order effects, and Soufi et al.|1998 for third-order 



use have shown their limits ( Reese et al. 2006 Suarez et 
. 20101. On the one hand these perturbative methods 
approximate the effects of the centrifugal force on acoustic 
pulsations in stars that show very high surface velocities, 
such as 5 Scuti and Be stars. On the other hand, they ap- 
proximate the impact of the Coriolis force on gravity modes 
in stars whose surface rotates slowly, but in which the pulsa- 
tion periods are of the same order as their rotation period 
(Prot ~ -Pose), such as SPB stars or 7 Doradus. Finally, 
perturbative methods may also fail in modeling pulsations 
of cooler stars, such as sub-giant or red giant stars, with 
very low surface velocities but rapidly rotating cores, and 
in which all pulsation modes are of a mixed nature, i.e. 
gravity close to the core and a coustic in the envelope (see 
for example 



Beck et al. 2012). 



effects). But these methods, although elegant and simple of 



This concerns many stars in C0R0T and Kepler obser- 
vation fields. Hence, if one wants to correctly extract the 
rotation profile from seismic observations, they need to cor- 
rectly apprehend the effects of rotation on pulsations. 

This work aims at presenting a model which accurately 
takes into account the non-perturbative effects of rota- 
tion on oscillation spectra. A new two dimensional non- 
perturbative code is presented. The 2D non-perturbative 
calculations fully take into account the centrifugal distor- 
tion of the star, and include the full influence of the Coriolis 
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acceleration. This 2D non-perturbative code is useful for 
studying pulsational spectra of highly distorted evolved 
models of stars, as well as stars presenting highly differ- 
ential rotation profiles. Section [2] introduces the basic pul- 
sation equations in spheroidal geometry using a coordinate 
system adapted to the star's shape. Section [3] is dedicated 
to the numerical method which is based on a spectral ap- 
proach for the angular part of the modes, and a finite differ- 
ence method, which is accurate to fourth order, for the ra- 
dial part. In Sect.[4j we test these calculations and evaluate 
their accuracy. Finally, in Sect.[5j t he results are vali dated 



by comparing them with those of |Reese et al. (20061, and 



an example of application is given in Sect. |6j 



2. Basic equations in spheroidal geometry 

When dealing with the subject of computing the pulsations 
of rotating stars, one has to face two main issues. Firstly, 
rotation, through the centrifugal force, distorts the reso- 
nant cavity of the pulsations. If a solenoidal rotation pro- 
file is assumed (around a north-south axis), the azimuthal 
symmetry is conserved whereas the spherical one is broken: 
the star takes on a spheroidal geometry. This centrifugal 
distortion, if it is strong enough has to be treated by a 
two-dimensional approach. Secondly, when rotation is ac- 
counted for in the pulsation equations, the Coriolis accel- 
eration enters the momentum equation and modifies the 
dynamics of the modes. If small enough, this Coriolis ef- 
fect can be approximated by perturbative methods. But for 
moderate to high rotational velocities, as well as for high or- 
der g modes when P rot ~ Pose, the perturbative treatment 
is no longer relevant, and a non perturbative approach has 
to be adopted. 



2.1. Spheroidal geometry 

Given the distorted shape of a rotating star, we choose a 
new coordinate system which defines the star's surface at 
a constant pseudo-radial coordinate. To do so, we adopt a 
multidomain approach, which consists in dividing the phys- 
ical space into domains whose boundaries correspond to the 



model's discontinuities (e.g. Canuto et al. 1988): one do- 
main Vi , which encloses the star, and one external domain, 
V2, bounded by the stellar surface and a sphere of twice 



the equatorial radius. Following |Bonazzola et al. (1998), 
we introduce a coordinate system which goes from spherical 
symmetry at the center to a spheroidal shape at the stellar 
surface, and back to a spherical geometry at the external 
boundary of V2. In this system, the radial coordinate, £, is 
no longer r, the distance to the center. However, for a fixed 
colatitude 9, r and £ are related thanks to the following 
continuous and bijective function: 
In domain V\: 



r(£0) = (l-e)C+ 2 (RM 



(1) 



where £ ranges from to 1, e = l — R po i/R eq is the flatness, 
and r(Q = 1,6) = Rs(8) the stellar surface. 
In domain V%: 

r((,8) = 2e+ (1 - e) C 

+ (2C 3 - 9C 2 + 12C - 4) (R s (6)-l + e) , (2) 




Fig. 1. Coordinate system used in ACOR. Vi extends from 
the center to the stellar surface, and V2 encompasses the 
star. 



where £ ranges from 1 to 2, £ = 2 corresponding to a spher- 
ical surface which encompasses the star (see Fig. [I]) . The 
use of such a coordinate system helps significantly with es- 
tablishing the boundary conditions. 

Then one has to define a new vector basis. We choose the 
following orthogonal spheroidal basis: 
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(3) 



where 



• a'e is tangential to iso-£ surfaces in the meridional plane, 



• a 



is directly orthogonal to <zg in the meridional plane, 



and where rg = dgr. 

In a forthcoming paper, this method will be applied to 
a realistic stellar model using three domains: one which 
encloses the convective core, the second the radiative enve- 
lope, and the third an external domain, which allows us to 
avoid discontinuities at the top of the convective zone. 

2.2. Basic equations 

Here, we consider the stellar interior to be an inviscid self- 
gravitating fluid, where the magnetic field is neglected. 
Therefore its physics is governed by the conservation of 
mass, momentum and energy, the energy transfer equation, 
and Poisson 's equation for the gravitation al potential (see 
for example Kippenhahn & Weigert 1994). 



Using the Eulerian formalism, we compute the oscilla- 
tion modes as the adiabatic response of the structure to 
small perturbations - i.e. of the densit y, pressure, gravita - 
tional potential and velocity field. As in Unno et al. ( 1989 ), 
we perturb the hydrodynamics equations around the equi- 
librium state. Considering that the velocity field in the equi- 
librium state is only due to rotation, the linearized equa- 
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tions in the inertial frame are given by 



dtp' + v- Oqv' + Ao) = o, (4) 

Pt + nd v ) v'il e ; + 2fl x v' + (v'.Vn) r sin 9e v = 



+ 00 



-— Vp' - V*'+ ^Vpo- 
Po Po 



(5) 



(d t + v • V) 



p 



_ v .(v*_vm = o,(6) 

po ^ipoJ V Po r^oy 

V 2 $' = AnGp' , (7) 



where primed quantities denote Eulerian perturbations, 
whereas quantities with the subscript "0" correspond to 
the stationary state. The symbol e; stands for the spheri- 
cal basis vectors. Note that the energy conservation equa- 
tion has been replaced by the adiabatic relation. Given that 
the equilibrium state is stationary and axisymmetric, the 
time and azimuthal dependences are of the form e J M+ m ¥>) ; 
where u) is the pulsation frequency, and m the azimuthal 
order of the pulsation mode. 

We put the system of equations in non-dimensional form 
using the following transformations: 



M 



GM 2 



( Req 
\GM 



fin 



n = 



ilk 



where Qk stands for the Keplerian critical velocity, i.e. the 
rotation velocity at which the centrifugal force compensates 
gravity at the equator. 

From now on, we work in terms of non-dimensional vari- 
ables and drop the tilded notation. 

Rather than projecting the motion equation onto the 
basis vectors, we choose to decompose it into one radial, one 
poloidal and one toroidal field. This decomposition allows 
us to obtain separate partial differential equations for the 
radial and angular coordinates, and helps to reduce the 
number of unknowns, as will be shown later. 

Moreover, we apply the change of variable 7r' = p' / p so 
as to avoid singularity problems at the surface of polytropic 
models, and we introduce an auxiliary variable d$', as well 
as the equation relating it to $': d$' = d^/dC, = 3^$', 
thereby reducing Poisson's equation from a second order 
differential equation to two first order ones. 

2.3. Spectral expansion 

Thanks to an appropriate choice of the coordinate sys- 
tem, the behavior of the eigenfunctions on iso-C surfaces is 
smooth. Therefore, the angular behavior of pulsation modes 
is well described in terms of an expansion on the basis of 
spherical harmonics. Such a spectral method is known to be 
very well suited in fluid dynamics. For instance, for a func- 
tion of class C°°, the numerical error decreases as e~ aAI , 
where M is the number o f spherical harmonic s in the ex- 



pansion and a a constant ( Canuto et al. , 1988 1 



Therefore, we perform a spectral expansion of the un- 



know ns in terms of the spherical harmonics Y™ ( Rieutord 
1987). Any vector field can be decomposed into a radial, a 
poloidal and a toroidal part. Therefore, the components of 
the pulsation velocity field are expressed as follows in the 
vector basis (a^,ae,e v ) 



v' z =i MOYt,m(P,<P), ( 8 ) 

t>\m\ 
+00 

v' e = i (Vi(C)d e Y( im (e,ip) + wt(C)-^-jYt tm (9,(p)) , 
t-^ V smfci / 

t>\m\ 

+00 

< = - E {MQ^Ye,m(0,<P) + MC)doYt, m (0,<pj) ■ 

t>\m\ 

All other scalar unknowns, namely and p' , are ex- 

panded in the same way as the scaled pressure perturbation: 

+00 

*'=J2 */(oiW0,vo (9) 

£>\m\ 

2.4. Symmetries and mode classification 

Because of the symmetries of the equilibrium model with 
respect to the rotation axis and to the equator, one ob- 
tains a separate eigenvalue problem for each value of the 
azimuthal order, m, and each parity, par, with respect to 
the equator. Thus, for a given value of m, there are two in- 
dependent sets of ODEs coupling the spectral coefficients, 
one with I of the same parity as m, and the other with 
opposite parities. That means, that for a given m, when 
including M terms in the spectral expansion, Vj G [1, M], 



I = I m I +2(J - 1) + par, for ir' e , cj>' e , d<j)' e , m, ve,p' e , 



tp = I m I +2(J — 1) + 1 - par for wo 



(10) 
(11) 



with par = if I is of the same parity as m (even mode) , 
and par = 1 otherwise (odd mode). We then obtain a sys- 
tem of ordinary differential equations (ODE) of the variable 
£ for the coefficients of the spherical harmonic expansion 
U£, v e , wi p , n^, f/ t , & e , and d$' e . 

2.5. Projections 

After having expanded the unknowns on the basis of spher- 
ical harmonics, the second step consists in projecting the 
equation system onto the spherical harmonics basis which 
is truncated to M terms, as is the spectral expansion. In 
this subsection, the unknowns are generically designated by 
Xi(C, 9, tp). Each partial differential equation of the form 

E(Xi(C,e,p)) = E f£x iA (C)Y 4 , m (^^)j - (12) 

is replaced by a system of M differential equations in f, 
obtained by projecting these equations on a basis of M 
spherical harmonics: 



sin 



-in 



E hTx i/2 (C)Y, 2 , m (^) Y} um = 0, 



J =1 



(13) 



in which l\ and I2 take on the values defined in Eqs.(10) 
or (THY 
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The equilibrium quantities, which are expressed as func- 
tions of £ and 9, are implicitly contained in the operator 
E. 

The impact of the basis dimension (M) on the precision 
of the computations will be discussed in detail in Sect. |4.l| 

2.6. Boundary conditions 

In order to complete the eigenvalue problem, it is necessary 
to specify a number of boundary conditions. The system of 
equations contains 4 sets of first order ODEs. Thus, we 
impose 4 boundary conditions. As the system is solved si- 
multaneously for all layers, two boundary conditions are 
imposed at the center, one at the stellar surface and one on 
the external spherical surface of V% (see Fig. [I]). 

Taking boundary conditions at the center is a delicate 
problem because of the coordinate system. It consists in 
imposing the regularity of the solutions at the center. To 
do so, we take the limits of the equations as C goes to zero. 
The different scalar unknowns behave as follows 

and the components of the velocity field obey 

u e = OiC 1 - 1 ) and u = 0(C), 
^ = 0(^-1) an d VQ = 0, 

This results in two algebraic relations between the un- 
knowns. A detailed explanation of these central boundary 
conditions is presented in Appendix [B| 

At the stellar surface, a stress free condition is imposed: 
5p = 0. The stellar surface is assumed to be on an isobar, 
at C = 1, thus the boundary condition corresponds to: 



i(m£! 



Po 



(14) 



where the subscript V stands for equilibrium quantities. 

The external condition for the gravitational potential 
consists in imposing that it vanishes at infinity. This can 
be achieved by matching the gravitational potential to a 
vacuum potential at £ = 2, i.e. on the spherical external 
surface V2. The advantage in using a spherical boundary is 
that the Laplace equation is separable for solutions of the 
form Xi(r,8,ip) — Xi(r)Y™(9,tp). This gives very simple 
conditions for a continuous match at the £ = 2 spherical 
boundary. Specifically, for £ > 2, this equation is decom- 
posed over the spherical harmonic basis, and each compo- 
nent solved separately. For each t, this yields two indepen- 
dent solutions: & e = Ar e which diverges at infinity and is 
therefore discarded, and &' e = B / r^ +1 ) which vanishes. 
Hence, the corresponding boundary condition is 



d$'» = 



r 



(15) 



where = d^r. 



3. Numerical method 



The sp e ctral f orm o f equa tions (A.ll, (A. 6), (A. 9), (A.12), 
(A. 13), ( A.15 1, and (A. 16) constitute a first order ordinary 



differential system of 7 x M equations, in terms of the co- 
ordinate £. Among them, 4xM are ordinary differential 
equations for t he spectral coeffi cients ug, ir' e , $J, and d&' e 
(i.e. Eqs. A.l A. 13 A.15 and A. 16), and the remaining 



Draic equations for the spectral 



3 x M equations are alge 
coefficients vg, wg p and p' e (i.e. Eqs. A. 6 A. 9 and A.12). 

We choose to solve this system by an Newton- like 
method, which consists in taking a guess at the pulsation 
frequency Co and looking for the pulsation mode with the 
closest frequency to this guess. Solving the system yields a 
deviation, 5a, from the initial guess, a = a a + 5a is taken 
as a new guess for the next iteration, and this process is it- 
erated till convergence (three steps are in general enough). 

We therefore isolate the terms proportional to 5a, and 
obtain the following system of equations which can be put 
under the form of a matrix: 



dyi 



= (A n + 5aA 12 ) Vl + (A 21 + 5aA 22 )z x , (16) 



= (B11 + 5aB 12 ) Vl + (B 21 + 5a B 22 ) z x , (17) 

where y\ and Z\ are column vectors containing the spectral 
coefficients of the unknowns 



and 



(18) 



and where I = \ m \ +2(j — 1) + par, Vj G [1, M],. 
An, Ai 2 , A 2 i and A 22 correspond to the following equa- 
tions 

pseudo-radial motion (Eq. 
definition of d<&' (Eq. 



Poisson (Eq. |A.i6|) 



A. 



A.l) 
W 



continuity (Eq. |A.13[ ) 



(19) 



whereas Bn, Bi 2 , B 2i and B 12 correspond to the algebraic 
equations 



poloidal motion (Eq. 
toroidal motion (Eq. 
adiabatic relation (Eq. 



A.6) 
A~9) 



|AT2| ) 



(20) 



Thanks to the three last equations, the non differentiated 
unknowns can be express ed in terms of the differentiated 
ones with the help of Eq. (17). To do so, the matrix (B 2 i + 
5a B 22 ), which is a real square matrix of rank 3M, has to 
be inverted. It is straightforward to show that if 5a is small 
enough, 



(B 2 



5aB 22 y x = B21 1 - 5a B^ 1 B 22 B^ 1 + o(5a 2 



( 21 ) 

Thus, the matrix system can be written under the form: 

{A + 5aA Sa ) Vl (22) 
5aB Sa ) Vl (23) 



dyi 
zi = (B 



3.1. Radial discretization 

In the radial direction, structural quantities, and as a conse- 
quence eigenmodes, undergo sharp variations (for instance 
at the top of a convective core, or at the star's surface). 
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Therefore, in the radial direction, spectral methods are in- 
appropriate when dealing with realistic stellar models. We 
choose to discretize the differential equations using a fi- 
nite differences approach. The continuous domain of inte- 
gration is replaced by a discrete set of N r radial points. The 
number of points in the radial grid is an important factor 
which affects the numerical precision, as discussed in de- 
tail in Sect. 14. l[ Another characteristic of the discretization 
scheme, which comes into play in the numerical precision is 
the estimation of the derivatives. With classical finite dif- 
ferences methods, the more precise you get the less stable 
are the computations. 

We have adopted a fourth-order difference scheme, 
which relies on the following identity 

h h 2 h h 2 

yi + 2 y'i + ^y" = - 2 y ' i+i + 12 + ^ ' (24) 

where the primes denote derivatives with respect to £, h is 
given by h = Ci+i ~Qi, and the subscripts and "i+i" de- 
note the layer indexes. The great advantage of this scheme 
is that it is accurate up to /i 4 , while retaining high nu- 
merical stability, as it only involves two consecutive grid 
points. This finite differences schem e has already been im- 
plemented by Scufiaire et al. (2008) in the Liege Oscillation 
Code which has been proven to be very stable and accurate 
for the calculations of all types of pulsation modes in all 
kinds of stars. 

From now on, we drop the subscripts "i" from y\. 
Thanks to Eq. (22), it is possible to express the derivatives 



of y. The identity (24) can then be valid as long as the 
matrices coefficients in A and A& a are continuous and have 
continuous derivatives. We then get the following matrix 
equation at each layer i 



a tv* - 


a i+l Vi+i = Sa [/3+y t + 


(3- +1 y t+ i]+o{h 5 ) , (25) 


where 






<4 = 


h h 2 , 2 
h + - 2 A l + -{A 2 + 




Pt = 


2 A Sa, i + Y2 0^* Ago, i 


+ Asa. iAi + A'g a i ^ , 


a; = 


h h 2 
Id - 2 A *+i + Y2 (^+1 + ' 


Pi = 


h A 

— 2 AScr, i+l 





+ Y2 ( Ai + l A Sa.i+l + Asa,i+lA l+1 + A' Sa l+1 ) , 

Here, Id is the square 4M x AM identity matrix. The system 
of equations over the whole stellar interior can then be put 
under the canonical form: 



Ay = SaAs^y, 

where the vector y has been introduced: 



y = 



yi 

!I2 



3JN. 



(layers 1, • • • N r ) 



(26) 



(27) 



We impose boundary conditions which can be expressed 
as algebraic relations, as explained in detail in Appendix 



We note that Eq. ( 26 1 is a generalization of the classical 



eigenvalue problem Ay = Xy. 
3.2. Inverse iteration algorithm 

In order to calculate the eigenvectors, y, and eigenvalues, 
da, we generalize the classical inverse iteration alg orithm 
to the mor e genera l probl em formulated in Eq. (26), as de- 
veloped by Dupret (2001 1 in the stellar pulsations context. 

The estimate at step k + 1 of the eigenvector, 34+i, is 
obtained from the estimate at step k and the initial guess, 
S (To j of the eigenvalue by the formula 

y k+1 = (A - 5a As^ 1 As* y k . (28) 

If we assume that A is inversible and that A^ 1 A& a is di- 
agonalizable, it is easy to prove that this algorithm has the 
same geometrical convergence rate as the classical inverse 
iteration algorithm. The inverse matrix is not explicitly cal- 
culated but we solve the system: 



(A - Sa Asa) 3Vm = Aga-yk- 



(29) 



To do so, we perform a LU factorization of the system with 
partial pivoting. Then, we iterate solving the 2 triangular 
systems at each step. Note that if the initial guess for the 
eigenvalue is good, the algorithm converges quickly towards 
the solution even with a bad eigenvector as an initial esti- 
mate. 

The corresponding eigenvalue is then determined by the 
generalization of the Rayleigh ratio 



8a 



y* Agq. Ay 
y* A^Ag^y' 



(30) 



where y* and A* &a are the Hermitian conjugates of y and 
Asa, res pectively. It can be easily shown that Sa 7 given by 
Eq. (30), minimizes: 



S 2 = y*(A* - 6aA* Sa )(A 
= \\(A- SaAsa) y || 



5aA S a)y 



(31) 



4. Tests and accuracy 



A and Asa- being block diagonal matrices. 



As mentioned in the introduction, in order to be integrated 
to a stellar modeling chain for massive computations, the 
program has been developed with a constant concern for 
simplicity and rapidity of use. In this section, we assess the 
role of numerical parameters in the convergence process to- 
ward an oscillation mode, and establish the computational 
performances with respect to these parameters. 
What takes up the most computing resources in ACOR is 
the integrations over 9 of the equation coefficients (Eq. 13). 
These coefficients are evaluated at each radial layers (for 
N r layers) by projecting the equations onto the spherical 
harmonics basis (of dimension M). Therefore, by assessing 
the role of the two parameters M and N r , it is possible to 
define the optimal values for a good compromise between 
computation time, memory resources and accuracy of the 
results. 

Note that there is no automatic method that allows us 
to identify the modes. In this work, we followed the progres- 
sion of the modes, as we gradually increased the rotation 



5 



R-M. Ouazzani, M-A. Dupret, D. R. Reese: Pulsations of rapidly rotating stars 




low order p mode 



10 15 20 25 30 35 40 

Spectral resolution M 

Fig. 2. Convergence as a function of the number of spheri- 
cal harmonics, taken as the relative frequency difference be- 
tween computations using consecutive spectral resolutions, 
M-l and M, for an n — 3 mode dominated by an I = 1, 
m = component, and for three different rotational veloc- 
ities: 18.5% Slfc in blue, 38% ilk in orange and 59% ilk in 
red. 



rate from zero to a high value. We used the kinetic energy 
distribution in the meridional plane, such as in Fig. [7j to 
correctly select the mode at each step. 

All the tests presented in this paper are made assuming 
uniformly rotating polytropes of polytropic index N = 3 
(polytropic exponent 7 = 4/3) with an adiabatic index of 

r 1 = 5/3. 

4.1. Convergence tests 

In the ideal case where equilibrium quantities would be 
C°°, the numerical errors would decrease as e~ aM . Rotation 
induces non-spherical profiles for the equilibrium quanti- 
ties, and causes the eigenfunctions to depart from a single 
spherical harmonic. Therefore, the higher the rotation rate, 
the stronger the deviations from sphericity, and the larger 
the spherical harmonic basis has to be, as illustrated in- 
deed, in Fig. [2] The convergence calculations illustrated in 
Fig. [2] show that convergence is reached for 7 terms for 
il = 18.5%il k , 16 terms for il = 37.9%il k and 25 terms for 
SI = 58.9%Jlfc in the spectral expansion. This figure also 
shows that the convergence reaches machine precision. 

Concerning the convergence with respect of the radial 
resolution, we present in Fig. [3] the worst (bottom) and the 
best case (top) . The higher the radial order n, the more nu- 
merous the nodes are in the eigenfunction and the higher 
the radial resolution has to be. These plots also show that 
using a non-regular radial grid, whose number of points 
increases as we go outwards, allows us to increase the ac- 
curacy of computations. This is due to the fact that the 
computed modes are acoustic modes, with consequently a 
small wavelength in the outer layers. 

4.2. Numerical resources 

Among all the operations performed by the code, calculat- 
ing the projection integrals (Eq.[l3]) is the most demanding 
in terms of computational time, whereas inverting the two 
matrices A and Ago- requires the most memory. In Table 
[l] are indicated the memory and time resources needed by 
the computations with respect to the parameters M and 
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Fig. 3. Convergence as a function of the radial resolution, 
taken as the relative frequency differences between com- 
putations using consecutive radial resolutions, r and (r+1) 
points, at three different rotational velocities: 18.5% in 
green, 37.9% fi^ in light blue and 58.9% ilk m dark blue. In 
both panels, the mode is dominated by the i = 2, m = 
component. The upper panel corresponds to an n — 3 mode 
and uses irregular grids, whereas the lower panel shows an 
71 = 9 mode calculated with even distributions of points. 



Table 1. Numerical resources - i.e. time in minutes and 
memory in MB or GB - used by ACOR with a spectral 
resolution M and a radial one of N r . 



N r 


M 


Matrix size 


Time 


memory 


2500 


40 


96 x 10 fa 


17 min 20 s 


6.4 GB 


2500 


20 


24 x 10 6 


2 min 30 s 


1.6 GB 


2500 


10 


6 x 10 6 


24 s 


500 MB 


1250 


40 


48 x 10 B 


8 min 


3.2 GB 


1250 


20 


12 x 10 6 


1 min 10 


850 MB 


1250 


10 


3 x 10 6 


12 s 


250 MB 


625 


40 


24 x 10 b 


3 min 50 


1.6 GB 


625 


20 


6 x 10 6 


33 s 


425 MB 


625 


10 


1.5 x 10 6 


6 s 


125 MB 



N r . Note that A and As a are bloc matrices, their size cor- 
responds to the number of non-zero elements they contain. 

5. Comparison with Reese et al. (2006) 

After the numerical tests presented in Sect. 4j the aim 
of this section is to validate ACOR's results by compar- 
ing them with those from the TOP code (Two-dimensional 
Oscillation Program). The TOP code has been developed 
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low frequency axi- symmetric modes 




axisymmetric modes 



high frequency axi-symmetric modes 



1x10" 




.g 1x10" 
CO 



1x10" 

0.2 0.3 

Fig. 4. Discrepancies between eigenfrequencies computed 
by ACOR and TOP as a function of the rotation angular 
velocity for five axisymmetric modes: Top: for two n = 
3 modes, corresponding to £ = 1 (light blue) and I = 2 
(orange) in the non rotating case and one n = 1 g mode 
with I = 2; Bottom: n = 9 modes, corresponding to £ = 1 
(dark blue) and I — 2 (red) . 



by Reese et al. ( 2006 ) for two dimensional polytropes as a 
first step. The approach is based on a two-dimensional spec- 
tral method which uses Tchebichev polynomials for the ra- 
dial dependence. We present here the comparison between 
TOP and ACOR results for identical polytropic models. 

Roughly, a polytrope is a simplified model for which 
one assumes an ad- hoc relation between density (po) and 
pressure (po)' 



Po 



(32) 



where K is the polytropic constant and 7 is called the poly- 
tropic exponent which has been taken as 4/3 here. The 
detailed method u sed to compute the ro tating polytropic 
models is given in Rieutord et al. (2005). This method is 
an iterative scheme based on a spectral decomposition us- 
ing Tchebichev polynomials for the radial dependence, and 
spherical harmonics for the horizontal one. We subsequently 
interpolate these models onto a radial grid which is appro- 
priate for finite differences. 

Concerning the parameters used for the calculations 
within the two codes, the angular spectral resolutions have 
been fixed so as to reach convergence. It therefore depends 
on the rotation velocity and varies from 10 terms in the 
harmonics expansion for low rotation rates to 25 for the 
highest ones. For TOP, the spectral resolution in terms of 
the Tchebichev polynomials varies from 50 to 80 terms. 
For ACOR, the radial resolution is fixed at N r =2000 grid 
points. 




Rotation rate D. I fl k 
non-axisymmetric modes 



a 



o 
D 
I 

E 




0.1 0.2 0.3 0.4 0.5 0.6 

Rotation rate Q. I fl k 

Fig. 5. Behavior of the eigenfrequencies computed by 
ACOR and TOP with respect of the rotational angular ve- 
locity for two multiplets; Top: centroid modes with domi- 
nant component I = 1, m = 0, and £ ~ 2, m ~ 0, Bottom: 
sectorial and tesseral modes, dominated by components 
I = 1, m = ±1 anc l £ = 2, m = ±1 and ±2. 



5.1. Eigenfunctions comparison 

The major effect of centrifugal distortion is the loss of 
spherical symmetry, which results in the coupling of spher- 
ical harmonics of different degrees to describe the horizon- 
tal behavior of a single mode. In Appendix [C] are given 
the contributions of dominant spherical harmonics in the 
spectral expansion of two modes: an odd mo de do minated 
by an £ = 1, m = component (see Fig. |C.1|) and an 



m — component (see Fig. 
even one dominated by t = 2, m = (see Fig. |C.2[). The 
solid lines correspond to the calculations done with ACOR 
and the dotted lines to those done with TOP. These plots 
clearly show that the evolution of the angular behavior of 
the modes with respect to rotation obtained by the two 
codes is very similar. This allows us to validate the analyt- 
ical calculations as well as the inverse iteration algorithm, 
which converges onto eigenfunctions, while eigenfrequencies 
are computed a posteriori through the minimization of Eq. 
(311. 



5.2. Eigenfrequencies comparison 

Concerning the comparisons of eigenfrequencies, Fig. [4] 
shows the frequency differences between the two codes, 
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Fig. 6. Avoided crossing, illustrated by plotting the fre- 
quency with respect to the rotation rate for two n = 2 
modes (referred to as mode #1 in red and mode #2 in 
blue), which, in the non-rotating case, are £ = 1 and £ = 5. 
The harmonic degrees given in the graph is the degree of 
the dominant component for each mode. 



This is illustrated in Fig. [6] by the evolution of the fre- 
quencies of two coupled modes with respect to the rotation 
rate. The two modes are of the same symmetry with m — 0, 
and £ = 1 or 5. As explained earlier, their frequencies can- 
not be degenerated, therefore the crossing is avoided, and as 
shown in Fig. [7] they progressively exchange angular char- 
acteristics. With the help of the distribution of the pres- 
sure perturbation in the meridional plane, we show that 
the mode with geometry dominated by £ = 1 at moderate 
rotation rates (mode #1), ends up with a dominant I = 5 
component at high rotation rates. 



Therefore, this work (see also Reese et al. , 2009 1 shows 
that modes in rapidly rotating stars can no longer be iden- 
tified only by one angular degree £. Actually, when rotation 
increases, the different components in the spectral expan- 
sion are more and more coupled by the non-spherical terms 
of the system of equations. As a consequence, each mode is 
composed of a mixture of spherical harmonics of the same 
symmetry, and it is not even possible to follow a mode con- 
sidering its dominant component, as it can change during 
an avoided crossing. 



for odd and even eigenmodes in the low (top) and high 
(bottom) frequency regimes. Globally, the discrepancies 
between results from the two codes is of the order of 
0.0001% - 0.08% (3 x 10~ 3 il k at most) for p modes, and 
0.5% for g modes (1 x 10~ 2 fife at most) which seems very 
satisfying considering that the two codes rely on different 
and independent computations, and in particular a different 
treatment for the central boundary conditions. Once more, 
this makes numerical programming mistakes unlikely in our 
calculations. 

One of the observational characteristics of rotation in 
seismology is the rotational generalized splitting, i.e. the 
frequency difference between two modes with the same 
radial order and harmonic degree, but with opposite az- 
imuthal orders. Figure [5] shows the evolution of the struc- 
ture of two multiplets with respect to the rotation velocity. 
The plots in the bottom panel show that the two codes 
find the same structure for the multiplets regardless of the 
rotation rate (from to 60% O^) and the symmetry class 
of modes. The impact of rotation on centroid modes (top 
panel) is also the same with the two methods. 

The agreement between the two oscillation programs 
developed separately, not only on the eigenfunctions and 
on frequencies, but also on the structure of the spectra, 
allows us to validate the approach adopted by ACOR. 



6. Illustration: avoided crossing 

In quantum mechanics, an avoided crossing (or level re- 
pulsion) occurs, for instance, in a two level system (|1) and 
1 2)), placed in a magnetic fi eld which acts differently on th e 



two levels (see for example Cohen- Tannoudji et al. 1973 ) 



When the two states are coupled, the levels repulse each 
other, since the system's energy cannot be degenerated. 

Accordingly, in a rotating star, pulsation mode frequen- 
cies tend to cross because the modes are not affected the 
sa me way by rotation (pa rticularly by the centrifugal force, 
see |Lignieres et al. 2006[ ) . As two modes of the same parity 



cannot have the same frequency, an avoided crossing occurs 
during which the two modes exchange their characteristics. 



7. Conclusion 

A new oscillation code, which computes non-radial adia- 
batic pulsations in rotating stars has been developed. The 
accuracy of the calculations has been achieved thanks to a 
hybrid method based on a spectral expansion on the spher- 
ical harmonics basis and a fourth-order finite differences 
scheme. The code has been tested and validated in the 
present study for polytropic models, but no assumptions 
were made in the implementation concerning the structure 
of the model used as input to ACOR. 

Although we've limited ourselves to barotropic models 
in this paper, it must be emphasized that our code is fully 
able to handle non-barotropic rapidly rotating models as 
will be presented in a forthcoming paper. This is an im- 
portant point for studying the pulsations of realistic stellar 
models such as those inclu ding rotation ally induced trans- 
port processes according to Zahn (1992). Indeed, such pro- 
cesses lead to shellular rotation profiles in radiative zones, 
which are as a consequence non barotropic. Moreover the 
radial treatment based on a finite differences method which 
is accurate up to the fourth order is particularly well suited 
for stellar models presenting sharp variations of the struc- 
tural quantities. 
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Appendix A: Equation system 

Rather than projecting the motion equation onto the basis vectors, we choose to decompose it into one radial, one poloidal 
and one toroidal field. We obtain one equation without derivatives relative to the angular coordinates by projecting the 
motion equation onto e r . Then by taking the divergence, and the pseudo-radial componant of the curl of the motion 
equation, we get two equations without radial derivatives. 

A.l. Radial motion 



dir' 

— = - i ( m O + a) Mfu£ - i (mO + a) sm6M^v' g (A.l) 
<9<f>' 

- smdM^v'^ - — + M p , p' + Mn> tt' , 
with, , 



M c 




M M 




M 


= -2Qr c , 


M p> 


1 8p\ 




p 2 tKJo„' 




1 d P \ 




P dC)e, v 



(A.2) 



where the colatitude 8 is replaced by p = cos 8, so rg = — s'mSr^. As stated in Sect. |2.2[ we applied the change of 
variable tt' = p' / p in order to avoid singularity problems at the surface of polytropic models. 

A.2. Poloidal motion 

Here we calculate the divergence of the horizontal component of the motion equation. To do so, one has to apply the 
operator rVj_-, where for any vector X = X^ + Xgcig + X^e^: 

Vx • X = — ~ (sin 8X g ) c + — £-{X v ). (A.3) 
r sin ft 88 ^ r sins' dtp 

Once applied to the motion equation, this gives 



(sin0 (mfl + <y)v' g ) r — — — - — (tlsmO ( — sm8 + cosdW^) 
^ sm 88 V r J 



sin 8 88 y v ' sin6» 88 V y r ' */ c sinfl 



2imflg ^cot# + — ^ v' g + img 



90 . 



= IrV + 1 f 2 $' + - — + - — 1 d _( sin6d P J\ + 1 d ( shx9d PA (A a) 

r r r 2 88 r 2 96* sin 8 88 \ rp 88 ) ( sin 8 88 \rp 2 88 ' ) f ' 1 ' j 

where g is given by g(£, /it) = 1/(1 + (1 — p 2 )r 2 /r 2 ), and the operator L 2 is the angular momentum operator: 



i 2 * = : — t: tt^ sin 



i 8 . i a 2 * 



sm 8 88 88 J f sin 2 9 dip 3 

(l-p 2 )^-) +2p^\ - j^-^r ~, (A.5) 



<9 2 *\ i s 2 4- 

+ ^ a^y*,- ~ (W^)^' 
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whose eigenf unctions are the spherical harmonics: L 2 YJ n (6, <p) — i (£+ 1) Y™{6, <j>). Thus, the equation becomes 



= (mfi + a) 



i d_ 
sinfl 80 



sm6 V ' p 



1 d 

D vfl ~r—j ^(sintfu ) 
sin 6* do v 



iD c v' c - ismOD^v'g - sinf-D^ - - L 2 ir' - 1 L 2 $' 



with, 



D„ 



<9$' 
dp 

= nt 

where 
-mg 



D pP ' + D PI1 



d£ 
dp 



dp 



i(C, M ) = 2( A1 -(l- / x)^) , 



20 ( 1 - — cot i 
r 



r c ( r 2 ' dC, r d6 



on 



_ am 



(1 + 5) 



ftgt 

1-p 2 



D P = 

D P , = 
D n = 

D^fj, — 



d ( 1 — fi 2 dp 



dfi \ rp 2 dfi 

(1 - M 2 ) dp 
rp 2 dp 
d ( ' 1 — p 2 dp 



dp \ rp dp 
1 - p 2 fldp 



r \p dp r 



(A.6) 



(A.7) 



A. 3. Toroidal motion 



The curl of the motion equation is (Unno et al. 1989) 
d 



dt 



n d 



ei + (uj' .V ft) r sin 6e,p + (v'.V) 2f2 - (212. V) v' + (W) 2f2 + V x ( -Vp 



0. 



(A. 



where we introduce the vorticity vector: u> = V x v. In order to find the toroidal component of the motion, we take the 



pseudo-radial component of the above equation (i.e. along ap 



= i(mfi + a) 



Id, „ . , 1 dv' 
— — (sin 01,;) - ^-h-o 1 
sm ay ^ sin ay 



ft 



dv 



+ z(mfi + o)—v' v + is\n9R v v' v 



C 
<9p 



-ft/, 



1 a 



i6> 50 



(smQv'g) + iR^ir' + iR p p' 



where. 



R ( = at + (l-p 2 )-^c + nt(d+f) 
i? M = - 1 b + c - 20 + fi(e + h) t 



R Cp = -n(l-p 2 ) (2+ T ^gl 



R WI = fttg 
mflt 

v = (1 - p 2 ) 



Rtt — — 



m dp 



rp dp 
m dp 
rp 2 dp 



(A.9) 
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t is given in Eq. (A.7) and a, b, c, d, e, / and h are defined by: 

r <9f2 , 2^ r u 
a = — — -(I-//) — JT 
r£ oC, r °p, 



dp 



d = 



, . t dr a 2 (1 — it 2 ) 9r u 
at, r o/i 



2- 

r 

(1-2.9) / 1 



r oQ r 



dr u 



h = 1 — g 
r 



\ r dp r 2 



A. 4. Adiabatic relation 

Let us define quantities analogous to the Schwarzschild discriminant in the two dimensional case: 



Ar = 



1 / <91np 



d\np 



Therefore, the adiabatic relation given in Eq. @ can be written under the following form: 



An = 

r. 



(A.10) 
(A.ll) 



i(mfi + cr) [ - p' - 7r' j - Aft>£ + sin A„ = , 
VP liP / 

where we have introduced: A^ = g A^ /r, and Aq = A^/r^ — (1 — ^i 2 ) r M A^/r. 



(A.12) 



A 5. Conservation of mass 

After linearization of Eq. Q, one obtains: 



ifmfi + cr) — + V • v' + v' • V In p = 
P 



With the help of Eq. ([6| , the equation can be rewritten: 



i(mQ + a) vr' + — v' • Vlnp + V.v' = 

?ip r a 



which finally gives: 



^ = -i(mn + a) r -^n> - g |-(sin^) - 

aC Tip rsinfl w; rsm9 dip 

dv 'c 

+ C c v f ( -&neC u i/„ + C c 



9/x 



with, 



r 



i p \ r z op, oQ J r 2r 2r M 



r \ 1 1 p op r 



(A.13) 
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A. 6. Poisson's equation 

Computing the Laplacien of <£>' in the new coordinate system, one obtains a dimensionless version of Poisson's equation: 



1 



d 2 <p> 1 d& 



re d fr g 



1 d fr 2 \ 



@C V r C / r 2 dC, \ r C/ 



'r 2 sinem V ~m) ~ r c r 2 dC,dO ' r 2 sin 2 9 d<t> 2 



2 • o 7^(— smt 
r z sin of) tq 



2r e d 2 $' 



1 



<9 2 $' 



(A.14) 



As mentioned in Sect. 2.2 we introduce a new variable and add an auxiliary equation in order to have a first order 
differential system: 

= 9 C $' (A. 15) 



with, 



P P > p' + P d v d*' 
9 ^ r C 



p, 



Pd*> = 



dpi 

(2 + 2/i^ 



P*'.L 2 (&) 



(A.16) 



r d/u, 



Pd*>n = 2(1 — p. ) -Hf- g 
-,9 



P P > = r 2 g 



(A.17) 



Eqs. ( A.l ), ( A.6I, (A.9|, (|A.12|), (|A.13|, (|A.15|) and (|A.16| are the 7 equations which make up the hrst-order differential 
equation system whicr 



r yields 2D pulsations of any rotating model 



In the second domain V2, i.e. in the vacuum, only Eqs. (A. 15 1 and (A.16) remain with p' = for Poisson's equation 



Appendix B: Central boundary conditions 



As explained in Sect. |2.6| we choose to impose two boundary conditions at the center. The goal of this section is to find 
two algebraic equations to replace differential ones at the center of the star. 

Table [BTT] shows the central behavior of the first terms in the spectral expansion of the unknowns. As can be seen in 
the table, the parity of the we coefficients is different from the parity of other spectral coefficients, 



Table B.l. Central behavior of the first terms in the spectral expansion in different cases: in blue are given the odd 
modes and in red the even modes. 
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The global parity of the mode is defined by the parity of ue, vg, ir' e , <E^, d$' e and p\. Note that in the case of an 
odd axisymmetric mode, wq has no meaning since the spherical harmonic lf=o, m =o is spherically symmetric and has no 
toroidal component. A similar argument applies to the v component in even axisymmetric modes. 

In order to emphasize the dominant terms in the equations near the center, we propose to express the spectral 
components under the following form: 
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4(0 = cV(0 
MO = C'-'MO 



p4(0 - CV/(C) 
«*(0 = C^(C) 



(B.l) 



where the tilded quantities have a constant behavior at the center. For the derivative of the gravitational potential, we 
obtain: d<Z>' e = C £+1 d% + £( e ~ 1 6 f t . 

Introducing this variable change into the system of equations, and taking the limit when £ goes to zero, we obtain: 



General case: 



Radial motion (Eq. A.l) 



(mil + a) (1 - e) v! t - 2mfi(l-e)^ - i$' t - £tt' 1 - 2 (£ - 1) Jffl (1 - e) w'^ <^ _1 



M p pf t + M n Tt'i + 2 (I + 2) J™ 1 (1 - e) w' t+1 - d& e 



C 



e+i 



Poloidal motion (Eq. A. 6 1 



= 



((mfi + cr) ^(i + 1) - 2toO) ^ - 2mO u' t - 2 fi + 1) (£ - 1) j; n 
- -£(t+l)ir' e - -£(£+l)&^ c^ 1 - 2ne(e + 2)j^ 1 w' e+1 c e+1 



Toroidal motion (Eq. A. 9 1 
= 20 



lp - 1) (£ lp + 1) J£ + [l Xp + 1) JT lp < p -l 



+ [((mO + a)h p {£ lp + 1) - 2toO) w' hp - 2Q£ lp J£ +1 ({£ lp + 2) v' tlp+1 + v! llp+l 
Adiabatic relation (Eq A.12[ ) 



Continuity equation (Eq. A. 13 1 

-i d«/ 



C 



5C 



= [(-(*-l) +C ( )u' e + £(£ + l)v' e ]C 



1 1 "o 



Poisson's equation (Eq. A. 16) 



= (l-efC'p'z ~ (2^ + 3)C £ ^ 



where J™ is given by: 



I 2 — m 2 
AP-1 ' 



0. 



if € > |m| 
if £ < Iml 



(B.2) 



(B.3) 
(B.4) 

(B.5) 
(B.6) 
(B.7) 

(B.8) 



and where e is the flatness, given in Sect. |2.1| 

We note, first of all, that the equations of radial motion, of continuity as well as Poisson's equation are singular at 
the center. In order to enforce the regularity of those equations, the boundary conditions that are imposed, are given by 
the fact that singular terms go to zero when £ — > 0. Therefore the radial motion and Poisson's equations lead to 



=(mn + a)(l-e)u' e - 2mfi(l-e)^ - £& t - £n' e - 2 (I - 1) J} n O (1 - e) w' t _ x 
=(1 -e) 2 ~p' g - (2£ + 3)d& e 



(B.9) 
(B.10) 



Case | to | = and par — and £ = 0: 
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This case is specific in the sense that, near the center, uq — O (£). Thus, the radial motion equation is not singular, 
whereas the continuity equation is: 



c 



du' 



(l-e)a 



Po 

riP 



+ (c< i) 



Therefore, the algebraic equation to impose near the center when 

= -(l-e)(T: 



Pa ~, 



to |=0, par = and £ = is: 
(C c - 1) u' 



(B.ll) 



(B.12) 



Concerning Poisson's equation, it is also singular in this specific case, and Eq. (B.10) is still relevant. 



It seems then that we obtained the algebraic boundary cond itions necessary to impose regular behavior o f the 
unkn owns near the center. However, the condition given in Eq. (B.9 ) is redundant with a linear combination of Eqs. (B.3 1 
and (B.4). Hence, the system is under-determined at the center. It is not possible to enforce the right behavior to the 



unknowns near the center using the tilded variables. We therefore have to impose central conditions on the non-tilded 
spectral components directly. Here again, different cases have to be explored: 



General case: 

The general case is composed of the following symmetry cases: 

— | m | = and par = 1 , 

— | m and par = 0, 

as well as all the components after the first one for: 

— | m |= and par = 0, i.e. £ = 2, 4, 6, • • • 

— | to and par = 1, i.e. £ =| m | +3, | m | +5, | to \ +7, ■ ■ ■ 

(B.13) 
(B.14) 

£ £ 

2mQ(l - e)v e - - -ir' e (B.15) 

(B.16) 



The regularity of the velocity held leads to = ii£ — £ vg 

= w e 

The regularity of the motion equation leads to = (1 — e)(mQ + a)ue — 

£ 

The regularity of Poisson's equation leads to = d$^ — - $' e 



Case | m |= 0, par = and £ = 0: 

In this specific case, W£-i has no meaning, we thus impose 

= v , 

On 

T Q + 3 U 

J- 1 -t 

= d& , 

where the condition over the radial velocity component is obtained from the continuity equation. 



= (l-eK<7 2 ^^,\ 



(B.17) 
(B.18) 
(B.19) 



Case | to 0, par = 1 and £ =| m \ +1: 



The continuity equation, the toroidal and radial motion equations, and Poisson's equation give respectively 

= M| m | +1 — £ v\ m \ + x , 

= ((mfi + cr)\m \ (| m | +1) - 2mf2) w\ m] , 

- 2fi | m | J,™ |+1 ((I to | +2)v\ ml+1 + u\ m]+1 ) , 

= (mil + a) (1 - e) u\ m]+1 - 2mO (1 - e) v\ m]+1 - **j m|+1 , 

" Zn\ m \+i- 2 \ m \ ■'M+i n ( 1 -eMm|> 
I to I +1 



o = d$: 



M+i 



c 



M+i ■ 



(B.20) 

(B.21) 

(B.22) 
(B.23) 
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The components W£ are neglected at the center in the remaining cases. This is an approximation, but after verifications, 
this method is the most efficient and numerically suitable in order to enforce regular behavior of the solutions at the 
center. 



Appendix C: Comparison of eigenfunctions 
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Fig. C.l. Radial parts of the eigenfunction n' e for a low frequency mode dominated by an £ = 1, m = component, for 
an iV = 3 polytropic model uniformly rotating at 18.5%, 46%, 63.5% and 83.7% of the Keplerian break-up velocity. 
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Fig. C.2. Radial parts of the eigenfunction n' e for a low frequency mode dominated by an £ = 2, m = component, for 
an iV = 3 polytropic model uniformly rotating at 3.7%, 18.5%, 46% and 64.5% of the Keplerian break-up velocity. 
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